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We present results on bulk thermodynamic quantities as well as net baryon number, strangeness 
and electric charge fluctuations in QCD at non-zero density and temperature obtained from lattice 
calculations with almost physical quark masses for two values of the lattice cut-off aT = 1/4 and 
1/6 . We show that with our improved p4fa3-action the cut-off effects are under control when 
using lattices with a temporal extent of 6 or larger and that the contribution to the equation of state, 
which is due to a finite chemical potential is small for ji^/T <1. Moreover, at vanishing chemical 
potential, i.e. under conditions almost realized at RHIC and the LHC, quartic fluctuations of 
net baryon number and strangeness are large in a narrow temperature interval characterizing the 
transition region from the low to high temperature phase. At non-zero baryon number density, 
strangeness fluctuations are enhanced and correlated to fluctuations of the net baryon number. If 
strangeness is furthermore forced to vanish, as it may be the case in systems created in heavy ion 
colhsions, strangeness fluctuations are significantly smaller than baryon number fluctuations. 
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1. Introduction 

Heavy-ion collision experiments at RHIC and LHC lead to thermalized dense matter at small 
but non-zero baryon density, or equivalently chemical potential. Therefore it is necessary to study 
the bulk thermodynamics of QCD at finite chemical potentials. In this work, we use the Taylor 
expansion method [|l|] to study the equation of state, the number density and fluctuations of various 
quantum numbers on the lattice. We study 2+1 flavor QCD with tree level Symanzik-improved 
gauge action and p4fat3 -improved staggered fermion action [||]. The simulations are carried out on 
16^ X 4 and 24^ x 6 lattices on a line of constant physics with almost physical quark masses; the 
pion mass is about 220 MeV and the strange quark mass is adjusted to its physical value. We have 
scanned a temperature range approximately from 170 MeV to 500 MeV. We are using the exact 
RHMC algorithm [Q] to update configurations. Details on our simulation parameters can be found 
in[§. 

2. Taylor expansions of thermodynamic quantities 

For a large homogeneous system, the pressure of QCD with u, d and s quarks can be expressed 
as ^ 

^ = ^lnZ{V,T,^,„^,,^,), (2.1) 

where the partition function Z is a function of the volume V, temperature T and chemical potentials 
of u, d, and s quarks. We have not considered other species of quarks whose masses are much 
heavier. Due to the sign problem, the difficulty of a direct lattice calculation at non zero chemical 
potentials arises. We perform a Taylor expansion in terms of the chemical potentials 

and compute the coefficients Cijk at zero chemical potentials. When the sum / + j + k is odd, the 
coefficient Cijk is given as expectation value of purely imaginary operators and therefore vanishes 
exactly. This reflects the invariance of the QCD partition function under change of particle and anti- 
particle. The leading term cqoo gives the pressure at vanishing baryon density and can be calculated 
via the integral method. Results for the parameter values considered here have been presented in 
In this work, we will concentrate on the part of the pressure 

Ap = p{Ji)-p{n=0), (2.3) 

that arises due to non-zero chemical potentials, where p. = (/Xi,,/Xd,/i.v)- For i + j + k > 0, the 
coefficients 

are derivatives of the partition function, and can be calculated on the lattice, where fi = n/T. These 
coefficients provide information about other thermal quantities as well. For example, the strange 
quark number density expands in chemical potentials as 

|i = £ + 1) c,y(^+i)/i;,/lj/i^, (2.5) 
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and similarly the light up and down quark numbers. We can further consider fluctuations in these 
quantities. 

Alternatively, one can introduce chemical potentials for the conserved quantities baryon num- 
ber B, electric charge Q and strangeness S, which are related to Hi, , , via 



12 11 11 



(2.6) 



and compute e.g. the baryon density as 



«B = ^ («« + nd + n,) . (2.7) 

Then we can study densities and fluctuations in B, Q and S. 

In the following, we will regard u and d quarks as degenerate and consider 2 + 1 flavor QCD. 
With the definition }Xq = llu = jJ-d for the light quarks, the coefficients are 



(2.8) 

Ai=0 



where the subscripts denote the order of the derivative and the superscripts indicate the correspond- 
ing flavors. If not specified, the default superscripts will be qs and will often be left out. It is evident 



from Eqs. ( |2.6| ) that choosing /x„ = pid is equivalent to a vanishing electric charge potential /ig = 0. 



Now we discuss how to evaluate these coefficients on the lattice. Inserting Eq. ( ^ ) into 



Eq. ( [2.8| ) and integrating out the fermion fields in the partition function yields the coefficients as 
expectation values of operators that contain derivatives of the determinant of the fermion matrix M. 
For example the formula for C20 reads 



where A^^ and Na are temporal and spacial extent of the lattice, «/ is the number of quark flavors 
in question (here n/ = 2), and (• • • ) indicates taking the thermal average over the ensemble. On 
each configuration, derivatives of IndetM need to be evaluated up to the same order as the order of 
the expansion coefficients. These derivatives lead to the appearances of the inverse fermion matrix 
M^^ inside traces 

5 IndetM / ,dM\ 

Tr M"^— , (2.10) 



dfX \ dpi 

^^indetM „ / .d'^M\ ( .dM ,dM 



To avoid full matrix inversions, we use the random noise method in estimating such traces. Suppose 
we have generated a set of N random noise vectors R^"^ ,a = I, . . . ,N, then the trace can be estimated 
as 



Tr^M 



- ^ ^-Y R^'')^M-^R^"\ (2.12) 
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where is some arbitrary matrix. For each vector 7?^"^ only the linear system MX = R'^"^ needs 
to be solved. It is still quite expensive to compute all necessary operators, since a large number of 
random vectors is needed in order to get a satisfactory accuracy. Also, higher order coefficients are 
more expensive, because more operators are needed. For the 4th order coefficients one has 



C40 



1 j «/ / 5^1ndetM\ 



where 



.(nf\2 j d^XnA&m dXnA&mX ^(nfV^ I ( ^^IndetM^V 
/nf\^ I ^^indetM /5lndetM\^\ /nf\^ I /^IndetMV 

(uf I ^^indetMV (nf\-^ / /5lndetM\^\\ 



^^IndetM 



Tr M 



3Tr M 



dM 



-M 



dpi djX^ 



+2Tr ( 



a'^lndetM 



Tr M 



^d^M 



M 



_^dM ^_^dM\ 
dpi dfj. J 



4Tr Im' 



dM 



-M 



d^M 



-JYvlM^^—rM^^—r +12Tr M"^— M"^— M 

( ,dM ,dM ,dM ,dM 
-6Tr M"^— M"^— M"^— M"^ — 
V on dji dji dji 



d^M 



(2.13) 



(2.14) 



(2.15) 



Five matrix inversions per random vector are necessary here, while for the 6th order, 12 matrix 
inversions are needed. Depending on quark mass, temperature and particular operator, different 
numbers of random vectors are needed to obtain that the errors arising from the stochastic estimator 
are smaller than or of the same magnitude as the statistical fluctuations within the ensemble. 



3. Pressure and densities 

In this section, we will first show results for the coefficients, then use them in computing 
pressure and quark number densities. 

In Fig. [D, we show the coefficients C200 and coo4 on both A^t = 4 and 6 lattices. C200, also known 
as the fluctuation in u (d) quark number density, increase rapidly through the phase transition 
region. As one can see, the lattice cut-off effect is small and seems to be under control. Results for 
C002 from = 8 lattices |Q] further support this statement. The fourth order coefficient C4 shows 
a pronounced peak around Tc. To compare the quark mass dependence, C200 and C002 for u and s 
quarks respectively are shown in Fig. ^. The slope is steeper for light than for the strange quarks, 
which indicates a stronger sensitivity to the chiral transition for lighter quark masses. We also show 
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Figure 1: C200 on the left and coo4on the right for = ^ and 6. The second order coefficients increase 
rapidly from confined phase to deconfined phase at around 200 MeV, while the fourth order ones develop 
a peak there. Stephen-Boltzmann Umits of the free case for the action that we use are marked for both 
quantities and matched very well in the high temperature region. 
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Figure 2: Second order coefficients C200 and C002 
for u and s quark respectively on A^^ = 6 lattice. 




150 200 250 300 350 400 450 
T/IVIeV 



Figure 3: cf^ on A^^ — 4 and 6 lattices. 



c\\ in Fig. |3[ which approaches zero from below in the high temperature limit. 

Combining all the measured coefficients, we obtain pressure and number density according to 
formula {12) and (2^). In Fig. 0, we show the pressure difference ^pjT^ and light quark number 
density Uq/T^ at finite light quark chemical potential but zero strange quark chemical potential 
jls = 0, up to the 4th order. This should be compared to the pressure at vanishing chemical potential 
[Q], which rises rapidly to a value of about p/T^ ?a 14 above the transition. The finite density 
contribution to the pressure adds to this less than 10% for piq/T < 1. 



4. Hadronic fluctuations at zero and non-zero chemical potential 

Fluctuations of charge densities hbs.q are related by the fluctuation dissipation theorem to the 
second derivatives of the partition function with respect to the corresponding chemical potentials 
Mb,5,2- Here B,S,Q denote baryon number, strangeness and electrical charge, respectively. Using 
Eqs. Q we can rearrange the expansion coefficients c"j[* of the pressure to get the coefficients of 
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Figure 4: Pressure Ap/T'^ and light quark number density n^/T^ at /ij = and jJ-q/T = 0.2 , 0.4, 0.6 and 
1.0. Small differences are observed between Nr ~ 4 and 6, especially when flg/T is small. Light quark 
number density seems to develop a peak around 200MeV when ^q/T increases. 
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Figure 5: Quadratic and quartic baryon number fluctuations at vanishing net density as function of temper- 
ature. Preliminary data from (2+l)-flavor simulations with almost realistic quark masses are compared with 
previous 2-flavor simulations [2]. Both results have been obtained on 16^ x 4 lattices. 



an expansion in I1b,s,q^ defined as 




(4.1) 



E.g., the following two relations hold for C2qq = cf and c^qq = 

cf = - (c|o + c^j + CI2) , cf = — (c^Q + + cfj + + cl^) . (4.2) 

In Fig. ^ we show the first two diagonal expansion coefficients in Hb/T as function of temper- 
ature, which can also be interpreted as the quadratic and quartic baryon number fluctuations. We 
compare our preliminary results for (2-1-1 )-flavor and almost realistic quark masses to earlier results 
with 2-flavor and a pion mass mj^ ~ lOOMeV [^. The normalization is such that in both cases the 
same Stefan-Boltzmann value for large temperatures is reached, i.e. we have divided by the num- 
ber of flavors. An obvious shift in the curves reflects the shift in the transition temperature from 
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Figure 6: Quadratic and quartic electric fluctuations at vanishing net density as function of temperature. Pre- 
Uminary data from (2+l)-flavor simulations with almost realistic quark masses are compared with previous 
2-flavor simulations [2]. Both results have been obtained on 16^ x 4 lattices. 



about 220 MeV to 200 MeV. Moreover the sudden change in the quadratic fluctuations is more pro- 
nounced for the smaller masses and the Stefan-Boltzmann value is reached faster. Correspondingly, 
the peak in the quartic fluctuations is higher for smaller masses. 

The expansion coefficients in jJ-s/T are identical to that in jJ-s/T - although the strangeness 
chemical potential differs from the strange quark chemical potential by a different sign - and are 
shown in Fig. [l] and |^ In Fig. ^ we show the first two diagonal expansion coefficients in pLq/T . 
The qualitative picture is very similar to PLb/T although the quark mass dependence of the peak 
height is significantly weaker. 

Using the expansion coefficients in IJ.b,s.q/T, one can construct hadronic fluctuations at non- 
zero baryon number density. Up to fourth order correction in Hb/T we have the following relations 
for baryon number, strangeness and electric charge fluctuations Xbsq^ 



T2 



2cf + 12cf f^f + - 



2ce + 24ef^y + 



T J 



f)' 



(4.3) 
(4.4) 
(4.5) 



In Fig. ^ we show baryon number and strangeness fluctuations at finite baryon number density. 
It is obvious that both quantities are developing a peak for increasing Hb/T. However, the peak 
in Xb is much more pronounced since this quantity eventually diverges at the critical point in the 
(r — /XB)-plane. As we anticipated from Fig. ^ the peak height in Xb is about twice as large as 
in earlier calculations with larger quark masses [Q]. Note that higher order corrections are still 
important, especially the position of the peak will be /Xg-dependent only by including the next 
higher order. This has to be analyzed in more detail and eventually will allow to limit the range of 
values for Hb/T where the leading order result is reliable. 



The off-diagonal coefficients in Eq. 4. 1 are usually connected to correlations between baryon 
number, strangeness and electrical charge. The correlation of baryon number and strangeness can 
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Figure 7: Baryon number and strangeness fluctuations at finite baryon number density, controlled by a finite 
baryon chemical potential. Results are correct up to fourth order corrections in chemical potential and have 
been obtained on 16^ x 4 lattices. 



0.45 




Figure 8: Correlation between baryon number and strangeness for several values of the baryon chemi- 
cal potential from 16^ x 4 lattices (left) and the linkage between baryon number and electric charge with 
strangeness respectively (right). On the right panel we compare or preliminary data (full symbols) to pre- 
viously obtained results from partially quenched calculations (open symbols) [7], both obtained on A^^ = 4 
lattices. 



be expressed in terms of expansion coefficients as 



^ {{nsns) - ins) {ns)) = cff + 3cff (^)' + , 



m 



(4.6) 



and is sliown in Fig. We find tliat also this quantity is developing a peak for increasing chemical 
potential, thus the enhanced correlations suggest the vicinity of a critical point. Another interesting 
quantity is the "linkage" of strangeness and baryon number or electric charge which is defined 



as Csx 



„sx 



/cj, where X = B,Q. It is known to be a robust quantity, i.e. the cut-off effects are 



small. In Fig. ^ (right) we compare our preliminary results with almost realistic quark masses with 
previously obtained partially quenched results and slightly larger light quark masses [^. The two 
calculations show good agreement, thus also the quenching and quark mass effects seem to be small 
in this quantity. Both results on correlation and linkage between the different quantum numbers 
suggest that the basic charges are carried by quasi-free quark directly above the transition. This 
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seems to rule out the existence of bound states as dominant degrees of freedom in this regime [||]. 

5. Conditions at heavy ion colliders and constrained densities 

In general, the pressure, or higher derivatives of the partition functions with respect to chemical 
potentials, are dependent on at least 3 variables Huxhs or equivalently Hb.s.q- So far we chose 
Hb > 0, while holding }Xs = }JLq = 0. To compare with experiment, for instance heavy ion collisions, 
the chemical potentials might need to be adjusted to meet the conditions of particular event-by- 
event fluctuation analyzes A very natural choice of the chemical potentials is to constrain the 



strange quark density to zero. Due to the existence of non zero off-diagonal coefficients in Eq. |4J 
we find an increasing strangeness with increasing jis, even for jis = 0. In heavy ion experiments the 
total strangeness is zero. Below we outline a procedure to constrain the net strange quark number 
density to zero, subsequently order by order in our /ig expansion. The procedure can be easily 
generalized to constrain other charge densities to arbitrary values. This might be of importance, 
since experiments are often restricted to certain rapidity windows, which may alter expectation 
values of charge densities. 

We can express the strange quark number density iris) in terms of the expansion coefficients 
of the pressure. Up to the 4th order, it reads 

ns = -ns{(lB,h) = -c\ftLB-2cl^(is-cf,iil-lcl^(il(is-^c'il^^^^^^ (5.1) 

where jX = n/T, which means that the strangeness chemical potential /I5 is no longer a free param- 
eter but depends on /ig, 

„BS \ ( 'y„BS„BS^ _ T,„BS„BS'^„BS , A^BS'^ „BS „BS _ A„BS^„BS \ 



Therefore, the formula for the pressure is modified to 



(5.2) 



An / ^BS2\ / „BS„B54 „BS^„BS „BS^„BS „BS„BS\ 

V <2) \ 16c^l 4c^l 2c^l y ^ ' 

which contains off-diagonal coefficients cii,ci3, etc. On the quark level those coefficients are 
generally small numbers since they are not present in the free theory. However, on the hadronic 
level they contain the diagonal strange quark coefficients which have - at least in leading order - a 
non-zero Stefan-Boltzmann limit. Hence the constraints ns = and jJ-s = lead to a quite different 
dependence of the pressure on jJ-s/T, as can be seen in Fig.^ (left). The difference is almost 
negligible, when performing an expansion in the light quark chemical potential Hq/T instead. It is 
interesting to mention that with the constraint tig = 0, the pressure expansion in jJ-q/T and Hb/T 
are identical up to a trivial factor between the two chemical potentials, i.e. the relation }JLb = ^i^q 
holds in this case and we have 

^p/T\^JT)l^^^^^p/T\^B|T,^Q = 0)l^^^. (5.4) 

We have also computed the constrained baryon number fluctuations at finite baryon chemical 
potential Xb- Qualitatively, the two cases of = and = are very similar. However, it is 
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Figure 9: The pressure Ap/T'^ up to the second order for both constraints as labeled (left) and the ratio 
■-A^Xb/Xb as explained in the text (right) for various values of Hb/T. The differences between the two 
constraints are of the order of 30% for both quantities. Results have been obtained on 24^ x 6 lattices (left) 
and 16^ x 4 lattices (right). 



interesting to remark that the two cases reach different Stefan-Boltzmann limits for high tempera- 
tures (r oo). Taking this into account we show in Fig. ^ (right) the ratio ^Xb/Xb> where ^ is 
the ratio of the corresponding Stefan-Boltzmann values. As one can see, the difference below Tc is 
as high as 30%. 

6. Summary and conclusions 

We have presented a method to rigorously compute corrections to bulk thermodynamic quan- 
tities at non vanishing chemical potential, by performing a Taylor expansion in n/T. Our new pre- 
liminary results improved previous calculations in many ways: we went to smaller quark masses, 
finer lattice spacings and 2-1-1 dynamical quark flavor. We also showed how to calculate various 
hadronic fluctuations, starting from a theory which naturally is formulated in terms of quark fields, 
as QCD is. The Taylor expansion method provides a variety of input to heavy ion phenomenology. 

Our findings are that the finite chemical potential contribution to the pressure is blow 10%, 
up to a chemical potential of Hb/T < 3 and that various hadronic fluctuations develop a peak with 
increasing baryon chemical potential. This seems to hold true also for strangeness fluctuations, 
although the peak is much less pronounced in this case. Correlations between strangeness and 
other charges increase as well when approaching the critical point. 
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